Free energy calculations of a proton transfer reaction by simulated tempering 
umbrella sampling first principles molecular dynamics simulations 
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A new simulated tempering method, which is referred to as simulated tempering umbrella sam- 
pling, for calculating the free energy of chemical reactions is proposed. First principles molec- 
ular dynamics simulations with this simulated tempering were performed in order to study the 
intramolecular proton transfer reaction of malonaldehyde in aqueous solution. Conformational sam- 
pling in reaction coordinate space can be easily enhanced with this method, and the free energy 
along a reaction coordinate can be calculated accurately. Moreover, the simulated tempering um- 
brella sampling provides trajectory data more efficiently than the conventional umbrella sampling 
method. 



Understanding chemical reactions by molecular simu- 
lations is a challenging problem, because a chemical re- 
action usually involves bond breaking, which cannot be 
treated by molecular simulations of classical mechanics 
based on force fields. We need to use the first principles 
(or ab initio) molecular dynamics (AIMD) methods to 
deal with bond breaking. However, the time span that 
can be studied by the first principles molecular simula- 
tions is very much limited, because their computational 
cost is much higher than that of the simulations with 
force fields. 

The molecular simulations based on quantum or clas- 
sical mechanics are hampered by the multiple-minimum 
problem, in which simulations tend to get trapped in 
the local-minimum states of free energy. One can use 
generalized-ensemble algorithms to overcome this diffi- 
culty (for a reviews, see, e.g., [H). Monte Carlo (MC) 
and molecular dynamics (MD) simulations based on 
generalized-ensemble algorithms have been widely per- 
formed for many molecular systems in order to have 
efficient conformational sampling. Three well-known 
generalized-ensemble algorithms are the multicanoni- 
cal algorithm (MUCA) @, |j| (for the MD version see 
Refs. 4, 5]), the replica-exchange method (REM) [G] (the 
method is also referred to as parallel tempering [7[ and 
for the MD version, which is referred to as REMD, see 
Ref. 8]), and simulated tempering (ST) UGH- 

Recently, general formulations for the multidimen- 
sional MUCA, REM, and ST have been given [Il|-[l3|. 
In this article we introduce a special realization of the 
generalized ST, which we refer to as Simulated Tem- 
pering Umbrella Sampling (STUS). This is a generaliza- 
tion of the Umbrella Sampling (US) method [14[ and is 



closely related to the Rep lica- Exchange Umbrella Sam- 
pling (REUS) in Ref. 0- 

Let us consider a system that consists of N atoms, 
where atom i has coordinate r,. We write the set of the 
coordinates of the atoms as r = {r\.r2, ■ ■ ■ ,r^}. The 
original potential energy function is represented by E(r), 
and the temperature by T. 

We propose a ST method in parameter space. The pa- 
rameter here stands for a label that specifies the poten- 
tial energy function. In this STUS method, M different 
restraint potential energy functions (or, umbrella poten- 
tials) VijVa,--' ,Vm are used. Introducing a reaction 
coordinate £, we define the umbrella potential function 
V m by 



V m (£(r)) = k m [£(r) - d r , 



(m = I,2,--. ,M), (I) 



where k m are force constants and d m are equilibrium dis- 
tances of the reaction coordinate. The STUS simulation 
yields a uniform probability distribution in parameter (or 
label) space. It means that a random walk in the M dif- 
ferent umbrella potential functions is realized during the 
simulation. 

In the STUS method, each state is specified by co- 
ordinate r and label to, and the following probability 
distribution function W m is used: 



W m (r) = Z- 1 cxp{-/? [E(r) + V m (r)] + a m } 
(to= 1,2,-- • ,M), 



(2) 



where /3 (= I/ZcbT 1 ) is the inverse temperature (&b is the 
Boltzmann constant), Z is defined by 



M 

E 

m— 1 



J dr cxp {-p [E(r) + V m {r)] + a m } . 



(3) 
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and introduced so that the probability distribu- 

tion in parameter space may be uniform. If the prob- 
ability in parameter space is constant, formally 
written as the following dimensionless free energy except 
for a constant: 



In jy dr cxp{~(3[E(r) 



V m (r)}} 



(4) 



We update the umbrella potential function to another 
one every few steps during the STUS simulation. When 
we attempt to change the m-th umbrella potential func- 
tion V m to the n-th umbrella potential function V n , we 
can use the following transition probability w so that the 
detailed balance condition may be satisfied: 



w(m — » n) — min ( 1, 



where 



Wn 



min (1, exp(— Z\)) , (5) 



A = j3[V n (r)-V m (r)]-(a ri 



0- 



(0) 



We have used the Metropolis criterion [16j] to satisfy the 
detailed balance condition. 

The ST parameters a rn (m = 1, 2, • • ■ , M) can be de- 
termined by the (multiple-histogram) reweighting tech- 
niques a ppli ed to a preliminary replica-exchange simula- 
tion |lll - ll3j . Once the ST parameters are determined, 
the STUS simulations are performed by repeating the 
following two steps: (1) perform a usual molecular sim- 
ulation with the potential energy function E + V rn for 
some steps, (2) update the umbrella potential V m to a 

I 



"neighboring" umbrella potential V n by the transition 
probability in Eqs. ([5]) and ([6]). If accepted, replace V m 
by V n . Go back to step (1). 

The free energy as a function of the reaction coordinate 
£, or the potential of mean force (PMF), J-(^o) for the 
original, unbiased system is defined by 



■Ffo) = -kvThxlZ- 1 / dr *(£(r) - exp[-/3£J(r)] 



(7) 

where Z is the partition function at temperature T, and 
8 is the delta function. Note that the umbrella potential 
functions are not included in this equation. 

We can use another expression for PMF instead of Eq. 



F(0 = ~k B T In P(0-C, 



(8) 



where P(£) is the probability distribution of the reaction 
coordinate £ in the original system without umbrella po- 
tential functions, and C is an arbitrary constant to set 
the zero point of the free energy. This form is conve- 
nient for molecular simulations because the probability 
P(£) can be obtained by generating a histogram of the 
reaction coordinate. 

From the results of the STUS simulation with the um- 
brella potential functions, we can calculate P(£) using 
some reweighting techniques such as the multistatc Ben- 
nett acceptance ratio (MBAR) estimator |17l|. which is 
based on the equations in Refs. [l8| and [l9j. 

The MBAR equations for calculating the expectation 
value (A) of a physical quantity A are written as follows: 



A(r n (k))exp[f - PE(r n {k))] 



M N n 

n ~ lk - 1 J2 N m exp{fm ~ P[E(r n (k)) + V m (r n (k))]} 



m—1 

M N n 



exp[-/3E{r n (k))} 



A I 



71=1 k = l 



£ N m exp{/ m - p[E(r n {k)) + V m (r n (k))]} 



m—1 



M Nn 



/* = -!*££ 



n=l k=l 



exp{-j?[£(r B (fc)) + V l {r n (k))}} 

M 

N m exp{/ m - (3[E{r n {k j) + V m {r n {k))]} 



-, (1 = 1,2,- •■ ,M), 



(9) 



(10) 



(11) 



where N m (m = 1, 2, • • • , M) are the total numbers of 
trajectory data in the simulation with V m , and r n (k) are 
the fc-th coordinate data in the trajectory obtained with 
V n . We can obtain the probability distribution P(£) of £ 
by calculating the expectation values of the number of £ 
taking the value £j (i = 1, 2, • • • ), where £j are discretized 



values along the reaction coordinate £. 

As an application of the STUS method, we considered 
the intramolecular proton transfer reaction of malonalde- 
hyde (see Fig. Q}. We performed STUS AIMD simula- 
tions of malonaldehyde in order to show that the STUS 
method is effective in calculating the free energy of chem- 
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FIG. 1. Intramolecular proton transfer reaction of malonaldehyde. The hydrogen atom that transfers between the two oxygen 
atoms is written as H. The two oxygen atoms that can bond the hydrogen atom are written as Oi and O2. 



ical reactions. 

We defined the reaction coordinate £ of the proton 
transfer as the difference of two distances between the 
hydrogen atom H and the two oxygen atoms Oi and O2: 



£M = |r Ql -t"h| - \ro 2 ~rn\ 



(12) 



When malonaldehyde is in a stable state where hydrogen 
atom H bonds to oxygen atom Oi, |ro x — ^h| should 
be less than \ro 2 — j*h| and therefore £ < 0. Likewise, 
when hydrogen atom H bonds to oxygen atom O2, £ > 0. 
When the reaction coordinate has a value which is nearly 
equal to 0, malonaldehyde should be in a transition state 
of the proton transfer reaction. 

We used 11 different umbrella potential functions in 
Eq. JTJ, that is, M = 11. k m were all set to be 0.01 
hartree-bohr -2 , and dijtfc, ••• ,<2ii were equally spaced 
between —1.0 A and 1.0 A, namely, they were set to be 
-1.0 A, -0.8 A, • • • , 1.0 A, respectively. 

We prepared the system of malonaldehyde in 71 wa- 
ter molecules with periodic boundary conditions. We 
used the CP2K program (version 2.1) [20| to perform 
AIMD simulations based on the density functional theory 
with the Born-Oppenheimer approximation. In the den- 
sity functional calculations, we used the Becke exchange 
functional 2l| and the Lee- Yang-Parr correlation func- 
tional [22|. The pseudo potential proposed by Goedecker, 
Teter, and Hutter [H S3 was used. A 280 Ry density 
grid was employed. We used the mixed Gaussian and 



plane waves approach 25j . We carried out the molecular 
simulations in the canonical ensemble. The simulation 
cell was set to be a cubic box (13.82 Ax 13.82 Ax 13.82 
A). The temperature of the simulation system was set to 
be 300 K. We used the Nose- Hoover chain method [26l — 
|28| as a constant temperature algorithm, set the number 
of chains to be 3, and set the time constant of the Nose- 
Hoover chain method to be 10 fs. The time step of these 
simulations was set to be 0.5 fs. The umbrella potential 
function was changed during the STUS simulation follow- 
ing the transition probability in Eqs. ([5]) and ([6]). In the 
STUS simulation we attempted to update the umbrella 
potential function to a neighboring one every 10 steps, 
that is, every 5 fs. We performed four independent simu- 
lations with the same conditions described above except 
for the initial velocities of the atoms. The data were 



stored each time just after the ST update attempts. Af- 
ter discarding thermalization steps, we obtained 214 ps 
simulation data altogether from the four STUS simula- 
tions. 

In order to show that the STUS method is more effi- 
cient than the conventional US method, we also carried 
out AIMD US simulations with a fixed umbrella poten- 
tial function using the same 11 umbrella potential func- 
tions as in the STUS simulations. Each of the 11 US 
simulations was performed for 20 ps and the data were 
stored every 5 fs. Actually, the ST parameters a m for 
the above STUS simulations were determined by apply- 
ing the MBAR reweighting techniques to the results of 
these 11 US simulations. 

The STUS method realizes a random walk in parame- 
ter space during the simulation. As a result, the reaction 
coordinate corresponding to the parameter can be sam- 
pled much more widely than in a conventional molecular 
simulation. Figure [2] shows the time series of the label of 
the umbrella potential functions and those of the reaction 
coordinate £ in one of the four STUS AIMD simulations. 
Through the simulation, the label of the umbrella func- 
tions largely fluctuated between 1 and 11 with almost 
equal probability. The reaction coordinate also fluctu- 
ated greatly, following the change of the umbrella poten- 
tial. Note that there is an expected strong correlation 
between the two graphs. The reaction coordinate took 
on values between around —1.5 A and 1.5 A. This means 
that malonaldehyde was able to experience not only the 
stable states, where hydrogen atom H bonds to either 
oxygen atom, but also the transition state. Essentially 
the same results were also obtained for the other three 
simulations. 

From the trajectory data of the STUS AIMD simula- 
tions, we calculated the PMF and the probability distri- 
bution of the reaction coordinate using the MBAR esti- 
mator in Eqs. ©, (jTOJ) , and (fTTj) . which are shown in 
Fig. [3] The PMF and probability distribution from the 
results of the conventional AIMD US simulations are also 
shown in Fig[3l 

One would expect that the free energy of the proton 
transfer reaction of malonaldehyde is a double- well func- 
tion which has minima at the two stable states and a lo- 
cal maximum at the transition state. Figure [3] (a) shows 



4 




FIG. 2. Time series of (a) the label of the umbrella potential functions and (b) the reaction coordinate £ during the STUS 
molecular dynamics simulations. 
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FIG. 3. (color online), (a) Potential of mean force (PMF) of the proton transfer reaction of malonaldehyde. The red filled 
circles represent the PMF obtained by the STUS simulations and the green filled squares by the conventional US simulations. 
The error bars were calculated by the MBAR estimator [l7| . (b) Probability distribution of the reaction coordinate £. The 
red curve is the probability obtained by the STUS simulations and the green curve is that by the conventional US simulations. 
The error bars are suppressed to aid the eye. 



that the PMF obtained from the STUS simulations has 
a simple double-well function which has two minima at 
nearly £ = —0.6 A and £ = 0.6 A, and a local maximum 
at around £ = 0.0 A. 

While the results of the STUS AIMD simulations were 
able to provide accurate PMF, the PMF calculated from 
the conventional US AIMD simulations could not be ob- 
tained accurately. We can understand this clearly com- 
paring the probability distribution of £ obtained by the 
STUS simulations to that by the US simulations. Al- 
though the probability in the two stable states is nearly 
equal to each other in the STUS simulations, the proba- 
bility obtained by the US simulations is not equally dis- 
tributed in the two stable states (see Fig. [3](b)). This 
implies that the US AIMD simulations could not sample 
sufficient trajectory data to calculate PMF in the time 
scale of the present simulations (20 ps for each simula- 
tion). 



In summary, we have proposed a new simulated tem- 
pering method, Simulated Tempering Umbrella Sampling 
(STUS), and applied it to the first principles molecular 
dynamics simulations of the intramolecular proton trans- 
fer reaction of malonaldehyde. We were able to obtain 
an accurate potential of mean force of the proton transfer 
reaction of malonaldehyde from the results of the STUS 
simulations. We also compared the potential of mean 
force obtained by the STUS simulations with the one ob- 
tained by the conventional US simulations. The STUS 
method is more efficient in exploring reaction coordinate 
space than the usual US method. 

In the present version of STUS, we fixed the tem- 
perature during the simulation. We can easily gener- 
alize STUS so that a two-dimensional random walk in 
both temperature and umbrella potential, as was done in 
REUS 

Moreover, the STUS method can be easily imple- 
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merited in the existing program packages. One does not 
need to modify the existing program. Because only the 
difference of energy is needed in the simulation, one has 
only to write a simple script that extracts values of energy 
and reaction coordinates and evaluates the ST transition 
probability during the simulation. 

Some of the computations were performed on the 
supercomputers at the Information Technology Center, 
Nagoya University, at the Research Center for Computa- 
tional Science, Institute for Molecular Science, and at the 
Supercomputer Center, Institute for Solid State Physics, 
University of Tokyo. This work was supported, in part, 
by Grants-in-Aid for JSPS Fellows, for Scientific Re- 
search on Innovative Areas ( "Fluctuations and Biological 
Functions"), and for the Computational Materials Sci- 
ence Initiative from the Ministry of Education, Culture, 
Sports, Science and Technology (MEXT), Japan. 
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